Variational Monte Carlo study of the ground state properties and 
vacancy formation energy of solid para-H2 using a shadow wave 

function 

Francesco Operetto 
Dipartimento di Fisica, Universita di Trento, 
via Sommarive, 14 1-38050 Povo, Trento, Italy 

Francesco Pederiva 
Dipartimento di Fisica, Universita di Trento, 
via Sommarive, 14 1-38050 Povo, Trento, Italy and 
INFM DEMOCRITOS National Simulation Center, Trieste, Italy 

(Dated: February 2, 2008) 

Abstract 

A Shadow Wave Function (SWF) is employed along with Variational Monte Carlo techniques 
to describe the ground state properties of solid molecular para-hydrogen. The study has been 
extended to densities below the equilibrium value, to obtain a parameterization of the SWF useful 
for the description of inhomogeneous phases. We also present an estimate of the vacancy formation 
energy as a function of the density, and discuss the importance of relaxation effects near the vacant 
site. 

PACS numbers: 67.80.Mg,05.30.Jp 
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I. INTRODUCTION 



The properties of molecular para-hydrogen (p-H2 ) have been the subject of intense 
experimental and theoretical investigation over the years. One main research direction is 
related to the existence of a metallic phase at very large pressuresQ, Q, Q, Q, Q, fl Q On 
the other hand, the Bosonic nature and the very small mass of P-H2 would naturally lead to 
the search for the existence of a superfluid phase. The strong H 2 -H 2 interaction, however, 
causes the zero-pressure ground state of p-H 2 to be a crystal at T=0. The only possibility of 
finding p-H 2 in a fluid phase at temperatures lower than 6K|8j], where a superfluid transition 
is believed to occur, comes from geometry-constrained or strongly inhomogeneous systems. 

It has been shown by Sinzingdre et al.j3| using Path-Integral Monte Carlo methods 
that small p-H 2 clusters with a number of molecules N<13 preserve a liquid character, and 

n 

show the signatures of superfluidity. Recently Levi and Mazzarello lO] proposed that the 
existence of such liquid clusters might depend on the difficulty for the system to go over the 
nucleation energy barrier in a time comparable to the lifetime of the clusters themselves. 
The determination of such energy barrier depends on the knowledge of a number of physical 
parameters, but in particular on the solid-liquid interface energy. It was recently shown that 
superfluid p-H 2 could also be found in two dimensional layers deposited on substrates such as 
graphite. Recent PIMC calculations showed that at full coverage the density of a monolayer 
is too large for allowing the system to have a stable liquid phase. However, the presence of 
alkali impurities, that exert a very weak attraction on the hydro gen molecules, can stabilize a 
disordered phase reducing the strong p-H 2 -p-H 2 interaction ll|, |l2| . The same effect can be 
studied reducing the coverage and looking at two-dimensional p-H 2 clusters at the surface. 
If the number of atoms in a single cluster is small enough (N < 30) the stable configuration 
is a puddle of liquid[13]. The occurrence of surface melting at the surface of a thick (5 to 7 
monolayers) layer of p-H 2 was also studied by PIMC techniques for temperatures > 3k}i^. 

Recent experiments performed by Grebenev et al.jl^J on a OCS-p-H 2 immersed in 4 He 
and mixed 4 He- 3 He clusters clearly show that the hydrogen coating the molecule undergoes 
a transition at a temperature between 0.38K and 0.15K. This transition can be inferred by 
the change in the measured rotational spectrum, and therefore of the momentum of inertia, 
when the OCS-p-H 2 complex is included in a pure 4 He cluster and when also 3 He is present. 
Recent simulation work confirms the occurrence of this transit ion [r|. 
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From the theoretical point of view the simulation approach is the most effective in study- 
ing inhomogeneous phases. In this paper we extend to the study of p-H 2 a technique that 
revealed extremely powerful in the description of inhomogeneous phases of 4 He with large 
numbers of particles: Variational Monte Carlo calculations based on Shadow Wave functions 
(SWF) jl^. fisj | . The main property of SWF is the capability of describing the crystalline and 
the disordered phase of a quantum system with the same translationally invariant form. It 
was therefore possible to study coexisting liquid-solid jl^ and liquid- vapor (^(J systems in 
4 He as well as finite systems like 4 He clusters with and without impurities[2l|. This method, 
although not exact, gives quantitatively sig nificant results, as for instance a realistic esti- 
mate of the interface energy of solid 4 He[l9|, in situations where other techniques, like Path 
Integral MC or Diffusion MC, are harder to apply. Some early calculations using SWF for 
small P-H2 clusters (N<7) were performed by Rama Krishna and Whaley[22]. However, 
their form of the wave function used cannot be immediately extended to bulk systems. In 
this paper we introduce a SWF optimized and parametrized in a form directly usable in 
calculations where the parameters become functions of the local density of the system 
With this wave function, we carried out a systematic study of the ground state properties 
of the bulk solid at densities up to 1.6 times the equilibrium density, and also to metastable 
bulk states below the equilibrium, beyond the limit of stability of the ordered phase. We 
also present results for an estimate of the vacancy formation energy in p-H 2 as a function of 
the crystal density. 

II. METHODS 
A. Hamiltonian 

Molecular para-hydrogen is here described as a collection of N point-like particles inter- 
acting via a two-body potential: 

# = -4 £ +£«fa)- (i) 



The model potential we chose is the Silvera- Goldman (SG)[23J potential, which takes the 
form: 

v{r) = A[v rep + (v att + v dd )f c (r)], (2) 
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r m (A) 3.41 
D 1.28 
a 10.923340 
b 10.098343 
7 0.4122340 

C 6 1.69550147 

C 8 0.71379389 

C 9 0.07468938 

Cio 0.38990868 

4 (K) 31.5763295 



TABLE I: Parameters of the SG potential 23]. All parameters are non-dimensional, except for r m , 
which is given in A, and A, which is given in K 

where, using f = r/r m . 



v a tt(r) 
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1.0 if r > £> r r 



The values of the parameters are given in table I. 

The Cg/f 9 term is introduced as an effective many-body force, assuming that the leading 
term is a triple-dipole contribution, which is usually introduced as a three-body Axilrod- 
Teller force. However, this parameterization is effective for studying isotropic properties of 
the system, and the replacement of the three-body term with a pair term is a second order 
effect. Venation. Monte Carlo captions by Nonnan at alfl 

, using a Bijl-Jastrow- 

Nosanow 2^ trial wave function, show that within this model the GS potential is in qual- 
itative agreement with experiments in a wide range of densities; other model interactions 
(*. the Scna^-Watts potentialQ or the Bne, at at. potentialft 

give better values 



around the equilibrium density, but become worse at higher pressure. 



B. Shadow Wave Function 



We use as trial solution for the N-molecules Hamiltonian at T=0K a Shadow Wave 
Function [17I. Iisj] . having the following general form: 



K(R, S)<j>(S) dS, 



(4) 



where R = {ri...rjv} and S = {si...Sjv} are the molecular and auxiliary ("shadow") 
degrees of freedom respectively The kernel K is written as the product of a Jastrow-pair 
wavefunction involving the molecular degrees of freedom, times a term coupling molecular 
and shadow positions: 



K(R, S)=Y[ exp 

i<j 



1 

2 



N 



J} exp -C(Ti - SiY 



(5) 



i=i 



The shadow degrees of freedom are also correlated by 0(s), a Jastrow product of the same 
form of the one used for the molecules. The correlation pseudopotential is the rescaled 
intermolecular interaction! 281: 



[s) = exp 



-SJ2 v ( as ij] 

i<j 



(6) 



The parameters b, C, 5, a appearing in the SWF are determined variationally, by minimizing 
the variance of the expectation value of the Hamiltonian 

(*(R)\(H - Eo) 2 *(R)) 



(AE 



(7) 



(*(R)MR)) ' 

As already mentioned, SWF provides a stable crystalline phase despite its translationally 
invariant form, due to the implicitly induced many-body correlations. The stable phase, 
ordered or disordered, is determined by the variational parameters minimizing the variance 
of the energy. 

For comparison we also performed VMC calculations with a standard Bijl-Jastrow- 



Nosanow wavefunction 
lattice positions: 



251 ]. in which the atoms are constrained to remain around given 



= exp 



v 



nexpf-Cfc-Ri 



(8) 



i=l 



5 



where the R are vectors of a lattice. This wave function has been also taken as importance 
function for the Diffusion Monte Carlo result that we computed at equilibrium density. 



C. Simulations 

All simulations have been performed for a cell of constant volume V containing a number 
N of molecules such that p = N/V, and imposing periodic boundary conditions in order 
to reduce finite-size effects. The simulations have been performed for two different crystal 
geometries: face-centered cubic (fee) and hexagonal close-packed (hep). Most of the results 
for the fee crystal were obtained using N=108 molecules arranged on 3x3x3 elementary cubic 
cells of side a ce u = (4/ p) 3). For simulations in the hep phase we considered N=180 molecules 
in a cell made up of 5x3x3 elementary cells. In order to check the magnitude of finite size 
effects (in particular for the vacancy formation energy) simulations in the fee crystal with 
N=256 (4x4x4 elementary cubic cells) and N=500 (5x5x5 elementary cubic cells) have also 
been performed. 

The molecule-molecule interaction is truncated at the edge of a sphere of radius equal 
to half of the side of the simulation box L. The contribution from the potential energy 
outside the sphere is estimated by assuming that the pair correlation function is constant 
for distances larger than L/2. 

The expectation values of the observables of interest have been computed by means of 
the Variational Monte Carlo method. When using SWF, the expectation value of a local 
operator O(R) is given by: 



where M is the normalization of the wave function. This integral can be computed sampling 
the joint probability distribution for the molecular and shadow degrees of freedom: 



As we illustrated in previous work[29j], the use of plain Metropolis sampling considerably 
slows down the convergence of the results, especially when the crystal includes a vacancy. 
This is due to the fact that the probability for a particle to have a given displacement is 
conditional on the position of the corresponding shadow degrees of freedom, to which it 





n(R, S, S') = K(R, S)<f>{S)K(R, S')<f>{S'). 



(10) 
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is connected by harmonic-like terms. The structure of the system is therefore analogous 
to that of a collection of trimers, whose single components can have only limited relative 
moves. By using pseudoforces (i.e. gradients of the function to be sampled) in sampling the 
probability distribution 7T, it is possible to move all the components of a trimer together. 
This was demonstrated to be a very powerful tool in the study of vacancies in solid 4 He, 
where it was shown that both relaxation of the density around the empty site and a mobility 



of the vacancy can be achieved 
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30|. The use of gradients also improves the convergence 



in the disordered phase, and in metastable crystal phases, where the low density might 
give rise to local disordering of the system. The optimization of the variational parameters 
appearing in the SWF was performed using the optimization program by C.J. Umrigar and P. 
Nightingale implementing a modified Levemberg-Marquardt algorithm. The optimization of 
the parameters of the shadow-shadow correlation pseudopotential of Eq. (6) is particularly 
demanding from the computational point of view. In fact the local energy does not depend, 
configuration by configuration, on the value of 5 and a, but on the statistical weight of 
the configuration n(R,S,S f ) only. The global computational cost of our SWF, including 
the optimization stage, is comparable with that of a DMC simulation for the same system. 
However, once a satisfactory parameterization has been obtained, successive simulations 
have a cost which is at least one order of magnitude less than DMC. For this reason SWF 
can be used for simulations including up to several thousands particles, as it was done for 
the study of the solid-liquid coexistence in 4 He. 



III. RESULTS 



A. Parameterization of SWF 



The optimized variational parameters in the SWF each have a different density depen- 
dence. The density dependence is exploited in the Local-Density dependent version of the 
SWF (LD-SWF)[l9j that is used for describing inhomogeneous systems, like the solid- vapor 
interface, with a single wavefunction. For such reason the optimization has to be carried out 
accurately in the region close to the equilibrium density (which is experimentally found at 
Po = 0.02595A~ 3 [2j), and lower. The parameters b, connected to the width of the correlation 
hole of the molecules, and C, giving the inverse of the mean square displacement between 
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a 1.2049±0.0872 
ai -21.961±2.984 
<5 -0.0047±0.0008 
61 0.3717±0.0370 

TABLE II: Values of the coefficients of Eq. (11) for the fits of the shadow-shadow pseudopotential 
parameters a and 5 as a function of the density. 

particles and shadows, show a weak dependence on the density. In particular the parame- 
ter C is nearly constant and equal to lA~ 2 throughout the density range considered. The 
parameter b tends to saturate at high densities around a value 3.35A, while at equilibrium 
density its value is b = 3.25A. For densities lower than the equilibrium density the parame- 
ter b shows some oscillations, contained within 5% of the value at equilibrium value. Such 
oscillations might be due to the fact that the system and autocorrelations between configu- 
rations become very large, making the optimization difficult. For such reasons a reasonable 
choice in the LD-SWF is to keep the parameters b and C constant at the value optimal 
at the equilibrium density. The remaining variational parameters a and 5 show a stronger 
density dependence. As already mentioned the energy at a given density has a much weaker 
dependence on the values of the shadow-shadow correlation parameters, due to the fact 
that <f>(S) does not enter in the estimate of the local energy. In Fig. 1 we show the results 
of the optimization. The errorbars indicate the range of values for which the variation in 
the computed variational energy is within two standard deviations from the optimal value. 
These parameters need to be made functions of the density in the LD-SWF. It is reasonable 
to assume a linear dependence of the parameters a and 5 on the density 

a(p) = a + aip 

5{p) = 5 + 5 lP . (11) 

The coefficients of the linear fits are given in Table II. The interpolations are also plotted in 
Fig. 1. The differences between energies computed with the optimal values of the parameters 
and the values for the LD-SWF parameterization less than 1%. 
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FIG. 1: Dependence of the SWF variational parameters a (upper panel) and 5 (lower panel) on 
the molecular density. Lines show the linear fits to the computed values. 
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E/N(K) 


SWF 


-86.679(3) 


JN 


-86.214(2) 


DMC 


-87.931(5) 


Exp. 


-93.5 



TABLE III: Ground state energies per particle at the density po = 0.02609A 3 obtained by VMC 
with SWF and JN, compared to the DMC result and the experimental value from Ref|2J 

B. Ground state properties 

The SWF has been used to compute the variational energies for molecular p-H 2 . The 
results are reported in tables III-IV. In table III we compare the binding energy for the 
density po = 0.02609A~3, (corresponding to a specific volume v = 23.08cm 3 mol -1 ) obtained 
by SWF, variational Monte Carlo with a Jastrow-Nosanow wave function (JN) and Diffusion 



24 As it 



Monte Carlo (DMC). We also compare to the experimental value reported in Ref. 
can be seen the energy obtained by SWF improve the Jastrow-Nosanow result. Nevertheless 
the difference between the SWF and DMC results is still larger than IK. The DMC result 
at this density shows a quite large discrepancy (about 6%) with respect to the available 
experimental results. This is in agreement with the findings of Norman et al., who point 
out that other model interactions give a better estimate of the binding energy near the 
equilibrium density. 

In Table IV we report the variational results obtained with our SWF in the disordered 
phase and in the crystalline phase, for which values for the fee and hep crystals are given. 
The results are also plotted in Fig. 2. The binding energy for the hep lattice is lower than 
for the fee lattice for densities p < 0.03170A -3 . 

As already mentioned, when using shadow wavef unctions, the stable phase (crystalline 
or disordered) is determined entirely from the variational principle. Therefore, given a set 
of variational parameters, we do not know a priori if the molecules will stay localized on 
a lattice. In order to determine the phase of the system, we need to compute a crystalline 
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TABLE IV: Total, potential, and kinetic energies per particle, in K, as a function of density from 
VMC-SWF calculations. For densities at which the system was found to be a crystal, the values 
are listed for the fee and hep lattices (obtained with N=108 and N=180 molecules respectively). 



p 


Eo/N 


TIN 

1 


V/N 








0050 


-17 400(16) 


4 139(15) 


-21 540(10) 








0.0080 


-29.348(13) 


8.202(16) 


-37.476(10) 








0.0110 


-41.390(9) 


14.269(14) 


-55.543(7) 








0.0150 


-55.929(5) 


24.037(5) 


-79.966(3) 








0.0172 


-63.800(7) 


34.019(9) 


-97.713(6) 












fee 






hep 






E /N 


TIN 

1 


V/N 


E /N 


T/N 


V/N 


0.01940 


-70.862(5) 


43.469(5) 


-114.332(7) 


-71.288(6) 


43.499(11) 


-114.787(11) 


0.02150 


-79.553(2) 


52.029(3) 


-131.582(3) 


-79.597(6) 


52.665(9) 


-132.262(7) 


0.02509 


-86.354(2) 


66.541(4) 


-152.895(4) 


-86.494(4) 


66.502(9) 


-152.996(8) 


0.02600 


-86.679(3) 


71.400(5) 


-158.079(5) 


-86.813(5) 


71.395(15) 


-158.209(15) 


0.02900 


-83.425(4) 


83.868(9) 


-167.293(11) 


-83.571(7) 


83.877(19) 


-167.448(21) 


0.03170 


-73.296(5) 


97.142(10) 


-170.438(12) 


-73.442(9) 


97.111(19) 


-170.553(23) 


0.03400 


-58.691(8) 


109.237(18) 


-167.927(18) 


-58.484(3) 


109.375(5) 


-167.858(5) 


0.03700 


-30.281(11) 122.924(22) 


-153.206(25) 


-30.139(2) 


122.976(4) 


-153.115(5) 


0.04000 


8.861(2) 


141.570(3) 


-132.708(4) 


8.931(2) 


141.450(5) 


-132.519(2) 



order parameter of this form: 

1 N M 

O = —Y Y |e- ik "' ri , (12) 

i=l a=l 

where the k a , a = 1...M are vectors of the reciprocal lattice for which the order is monitored. 
This quantity is exactly 1 if all the molecules are sitting on top of a lattice site, while in 
a disordered system its value is about 1 / y/N, where N is the number of particles. The 
value we find at the equilibrium density is (O) ~ 0.72 for the molecular degrees of freedom. 
This value slightly increases with the density of the system as the molecules become more 
and more localized. At the highest density considered we found (O) ~ 0.76. The limit of 
stability for the crystal phase is found at a density 0.0172A~ 3 < p < 0.0194A" 3 . 
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FIG. 2: Estimated energy per molecule in solid and liquid P-H2 • Filled circles: SWF results for the 
fee crystal; empty triangles: SWF results for the hep crystal; stars: SWF results in the disordered 
phase. The inset gives an expanded view around the equilibrium density. Filled triangle: JN result; 
diamond: DMC result. The dashed line is a plot of the fit with Eq. (13). The arrow indicates the 
equilibrium density obtained by the fit. 



In order to compute the value of the equilibrium density and pressure, it is convenient to 
fit the results for the ground state energy per particle by means of the following expression: 



E(p) = E + ap + bp^ 



(13) 



The coefficients for the solid phase in the fee and hep crystals are reported in table V. 
The estimate of the equilibrium density is in good agreement with the experimental finding 



exp 



p r = 0.02595A 3 2]. The pressure as a function of density is then obtained from the 



following expression: 



P{P) 



1 dE(p) 
p 2 dp 



(14) 



In Fig. 3 the computed pressure for the fee and hep solid phases is compared to the 
experimental curve extrapolated to T=0 3]. The SWF result is lower than the experiment 
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FIG. 3: Pressure as a function of density in solid P-H2 • Dotted line: SWF result computed from 
Eq. (15) from the fit of the ground state energy. Filled circles are experimental data from ref. [jj. 
Points are values computed in ref. 12411 with different intermolecular potentials (MS: Meyer and 
Schaefer 31|; SW: Schaefer and Watts 20 : SG: Silvera and Goldman^]; BEL: Buck et aZ.Q)- In 



the inset we expand the same data near the equilibrium density. 

near the equilibrium density, while the pressure tends to be overestimated at high densities. 
We also compare our results with results of Norman et al. [24| for the same Silvera- Goldman 
potential and for other intermolecular interactions. 
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FIG. 4: Pair distribution functions computed with SWF in molecular P-H2 • Solid lines: pdf of 
particles; dotted line: pdf for shadows, (a) p = 0.00800A~ 3 (disordered phase); (b) p = 0.02609A -3 
(fee crystal); (c) p = 0.02609A~ 3 (hep crystal); (d) p = 0.04000A- 3 (fee crystal). 



In Fig. 4 we plot the pair distribution function: 

g (r) = (£6(1^ -Tj\,r)), 

i+3 



(15) 



were () indicates the expectation value over the SWF. The same quantity can be defined for 
the shadow degrees of freedom: 



9 sir) = Q26(\si-Sj\,r)). 



(16) 



g{r) and g s (r) are displayed for a density below p , at p for the fee and hep crystals and at 
an higher density for the fee crystal. At the lowest density g(r) has the typical behavior that 
can be observed in a fluid, with a peak at the first shell of about 1.1. At larger densities, 
where the system is crystallized, it is possible to observe how the first peak increases to 
values around 2. The peaks corresponding to the second and third shell are nearly merging 
into each other. This feature, common to other quantum solids, is due to the wide dispersion 
of the molecules around the lattice sites. On the other hand, the distribution of shadows 
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presents, as expected, a much better defined structure. The peaks typical of the fee and 
hep structures, already well visible at density po become separated at higher densities. The 
shadow-shadow correlation induces a much stronger localization, and at higher densities the 
distribution becomes similar to that of a classical solid. 



C. Vacancies 



The vacancy formation energy at constant pressure in a system of N particles and iVj 
lattice sites at a given density p can be defined asj^, 

N - 1 

AE vac = E(V, N = N-l,N t = N)- ~^-E(V, N, N l = N). (17) 

The volumes V and V are related by the requirement that the density of the system remains 
constant, i.e., N/V = (N — 1)/V. In a quantum crystal AE vac has contributions from three 
different effects. The main contribution comes from the missing potential and kinetic energy 
due to the missing particle. At constant pressure and in the thermodynamic limit, this energy 
is a function of the pressure and of the potential energy per particle: 

A< c = -P{p) ~ ^f. (18) 

Another important contribution comes from the relaxation of the crystal around the empty 
site. This term lowers the vacancy formation energy due to the fact that the kinetic energy 
of the molecules surrounding the vacancy is reduced. This contribution is more important 
at lower densities, and can be well described only within a model that allows lattice sites to 
be displaced from their original positions, as happens for SWF where no a-priori lattice is 
assumed. The third contribution comes from the motion of the vacancy through the crystal. 
The bandwidth of the vacancy has recently been recently estimated to be 6 to 10K in 4 He 
by Galli and Reattoj^J. We expect this contribution to be smaller in p-H 2 , due to the 
higher degree of localization of the molecules. 

In Fig. 5 we report the results obtained for AE vac as a function of the density for the 
fee crystal (obtained with N=108 molecules), and for the hep crystal (obtained with N=180 
molecules). The values are compared with the corresponding estimate of AE® ac . As can be 
seen there is a strong dependence of AE vac on the density. For low densities the contribution 
from lattice relaxation, which can be estimated by the difference E re \ = AE™ — AE vac , is 
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FIG. 5: Vacancy formation energy AE vac from SWF calculations as a function of the molecular 
density. Circles: fee crystal (108/107 particles); squares: hep crystal (180/179 particles). The 
dotted line represents the estimated energy for a vacancy without relaxation effects, obtained 
from the fitted equation of state. In the inset we plot the vacancy formation energy at density 
p = 0.02609A" 3 as a function of N, the number of particles used in the simulation. 

E re i ~ 9K. This value decreases when the density increases. For high densities E rei does not 
show a systematic behavior. The estimate of the energy of a static vacancy depends on the 
potential energy per particle, which has a larger uncertainty at higher densities due to larger 
autocorrelations and stronger finite-size effects. Finite size effects also make the estimate 
of AE vac from the variational results harder to obtain. This can also be inferred from the 
fact that the results for the fee and hep lattices tend to depart from each other in a non 
systematic way. Their difference can be considered as a measure of the actual uncertainty 
on the computed value of AE vac The use of periodic boundary conditions implies that the 
presence of one vacancy in the simulation cell corresponds to a density 1/N of vacancies in 
the infinite crystal. The local deformation of the crystal due to the relaxation of the atoms 
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N AE vac (K) 



107 149(1) 
179 150(2) 
255 154(2) 
499 150(3) 

TABLE VI: Vacancy formation energy for p-H 2 at density p = 0.02600A~ 3 as a function of the 
number of particles used in the simulation. 

around the empty site might induce deformations to the next shells, which might translate 
into a vacancy- vacancy interaction. In order to check if such interaction affects the estimate 
of the vacancy formation energy, we performed simulations at different numbers of molecules 
at the equilibrium density (where relaxation effects are large), and recomputed the vacancy 
formation energy. Results are reported in table VI. As can be seen, the differences are 
within two standard deviations, suggesting that no sizeable effects of the relaxation can be 
observed beyond the second shell of neighbors (which is contained in the simulation box 
with 108 molecules). This is also confirmed by the fact that at density p no significant 
difference is observed for the vacancy formation energy in the fee and the hep crystals, that 
have the same structure for the first shell of neighbors. 

IV. CONCLUSIONS 

A shadow wavefunction has been devised to study ground state properties of molecular 
para-hydrogen, using the Silvera-Goldman interaction. The equation of state of the crys- 
talline fee and hep phases has been computed optimizing the variational parameters of the 
SWF, and fitting them as functions of the density. The estimated energy per particle is 
on average 10% higher than the available experimental data. However the estimate of the 
pressure and of the equilibrium density suggest that this SWF gives a realistic description 
of the system, allowing for future developments, in particular for as regards the study of 
the solid-vapor interface, and the study of clusters. The vacancy formation energy has been 
estimated in a wide range of densities. At densities around the equilibrium value relaxation 
effects appear to be about 5% of the total vacancy formation energy. 
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